#pragma once
#include <algorithm>
#include <array>
#include <cmath>
#include <eigen3/Eigen/Core>
#include <functional>
#include <iostream>
#include <vector>
using namespace std;

template<unsigned int dim, unsigned int deg>
class Spline
{
protected:
	// info
	//Due to compatibility with STL containers, the subscripts of all our variables and functions start with 0,which is different from ZJU's handouts.

	int m_N; // the number of knots
	vector<pair<double, array<double, dim>>> m_points;
	unsigned int m_dim = dim; // dimension == 1 if y=f(x); dimension == d if we fitting scatter points in d-dimension space 
	unsigned int m_deg = deg; //degree; default 3; S_(m_deg-1, m_deg) splain
	array<function<double(double)>, dim> m_s = {}; // the s(f;x) spline function
public:
	Spline() //default; not recommanded
	{
		m_N = 0;
		m_points = {};
	}
	Spline(vector<double> t, function<array<double, dim>(double)> f) //pass parameters based on independent variables and function expressions
	{
		m_N = t.size();
		vector<pair<double, array<double, dim>>> points;
		sort(t.begin(), t.end());
		for (auto it = t.begin(); it != t.end(); it++)
		{
			pair<double, array<double, dim>> temp;
			temp.first = *it;
			temp.second = f(*it);
			points.push_back(temp);
		}
		m_points = points;

		/*for (auto element : m_t)
			cout << element << " ";*/
	}
	Spline(vector<pair<double, array<double, dim>>> points) //directly transfer parameters based on the coordinates of points
	{
		m_N = points.size();
		sort(points.begin(), points.end());
		m_points = points;
	}

	void printpointsData()
	{
		for (auto it = m_points.begin(); it != m_points.end(); it++)
		{
			cout << "Argument: " << (*it).first << "; ";
			cout << "Point coordinates: (";
			for (int i = 0; i < dim; i++)
			{
				cout << (*it).second[i];
				if (i != dim - 1) cout << ", ";
			}
			cout << ")";
			cout << endl;
		}
	}
	void printSplinefunData(unsigned int prec = 100) //precision is the number of intercells in [a, b]
	{
		double a = m_points[0].first, b = m_points[m_N - 1].first;
		for (int i = 0; i <= prec; i++)
		{
			double t = a + i * ((b - a) / prec);

			cout /* << "Argument: " */<< t << "; ";
			//cout << "Point coordinates: (";
			for (int j = 0; j < dim; j++)
			{
				cout << m_s[j](t);
				if (j != dim - 1) cout << ", ";

			}
			//cout << ")";
			cout << endl;
		}
	}




	
};

template <unsigned int dim, unsigned int deg>
class ppSpline : public Spline<dim, deg>
{
private:
	//two additional constraints for S_(2, 3)
	/*
	* complete cubic spline : s'(f;a) = f'(a) and s'(f;b) = f'(b)
	* cubic spline with specified second derivatives at its end points: s"(f;a) = f"(a) and s"(f;b) = f"(b)
	* natural cubic spline 
	*/
	array<double, dim> y_a, y_b;
public:

	using Spline<dim, deg>::Spline;

	// the d-th dimension (i = 0, 1, ..., dim-1); the j-th point (j = 0, 1, ..., m_N-1)
	/*
	* The original intention was to use template parameters for d and j, but it is not allowed in the bind() syntax in <functional>
	*/
	double linearppFormpiece(double t, int d, int j) // n=1, t is in [t_j, t_j+1] (0<=j<=m_N-2)
	{
		if (j == this->m_N - 1)
		{
			return this->m_points[j].second[d];
		}
		else if (this->m_points[j].first == this->m_points[j + 1].first) // If the denominator is 0, then the entire score value is considered 0.
		{
			return this->m_points[j].second[d];
		}
		else // j>=0 && j<=this->m_N - 2 && this->m_points[j].first != this->m_points[j + 1].first
		{
			return (t - this->m_points[j].first) * (this->m_points[j + 1].second[d] - this->m_points[j].second[d]) / \
				(this->m_points[j + 1].first - this->m_points[j].first) + this->m_points[j].second[d];
		}
	}
	double linearppForm(double t, int d)
	{
		for (int j = 0; j < this->m_N; j++)
		{
			if (j != this->m_N - 1)
			{
				if (t >= this->m_points[j].first && t < this->m_points[j + 1].first)
				{
					return linearppFormpiece(t, d, j);
				}
			}
			else // j == this->m_N - 1, i.e. t is b
			{
				return linearppFormpiece(t, d, j);
			}
			
		}
	}


	//Only support cases where n=1 or 3
	//t is independent variable
	void ppForm()
	{
		if (deg != 1 && deg != 3)
		{
			cout << "illegal parameter!" << endl;
			exit(-114514);
		}
		else if (deg == 1) //linear interpolation
		{

			for (int d = 0; d < dim; d++)
			{
				this->m_s[d] = bind(&ppSpline::linearppForm, this, placeholders::_1, d);
			}
		}
		else // cubic interpolation
		{
			//To do...
		}

	}

};

